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Abstract 

A computer program is described to calculate momentum distributions in stripping and diffrac- 
tion dissociation reactions. A Glauber model is used with the scattering wavefunctions calculated 
in the eikonal approximation. The program is appropriate for knockout reactions at intermediate 
energy collisions ( 30 MeV < E^afe/nucleon < 2000 MeV). It is particularly useful for reactions 
involving unstable nuclear beams, or exotic nuclei (e.g. neutron-rich nuclei), and studies of single- 
particle occupancy probabilities (spectroscopic factors) and other related physical observables. 
Such studies are an essential part of the scientific program of radioactive beam facilities, as in for 
instance the proposed RIA (Rare Isotope Accelerator) facility in the US. ^ 

PACS numbers: 25.60.Gc, 24.50.+g, 25.60.-t 
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^ "Neutron saturated nuclei are the closest one can get to having a neutron star in the laboratory. The 
study of drip-line nuclei has progressed remarkably by observing nuclear reactions caused by radioactive 
fragments." P. Gregers Hansen, Nature 334, 194 (1988). 
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PROGRAM SUMMARY 



1. Title of program: MOMDIS (MOMentum Distributions) 

Computers: The code has been created on an IBM-PC, but also runs on UNIX or 
LINUX machines. 

Operating systems: WINDOWS or UNIX 

Program language used: Fortran- 77 

Memory required to execute with typical data: 16 Mbytes of RAM memory and 2 MB 
of hard disk space 

No. of lines in distributed program, including test data, etc.: 4320 
Distribution format: ASCII 

Keywords: Momentum distributions; Breakup, Stripping, Diffraction dissociation 

Nature of physical problem: The program calculates bound wavef unctions, eikonal 
S-matrices, total cross sections and momentum distributions of interest in nuclear 
knockout reactions at intermediate energies. 

Method of solution: Solves the radial Schrodinger equation for bound states. A Nu- 
merov integration is used outwardly and inwardly and a matching at the nuclear surface 
is done to obtain the energy and the bound state wavefunction with good accuracy. 
The S-matrices are obtained using eikonal wavefunctions and the "t-pp" method to ob- 
tain the eikonal phase-shifts. The momentum distributions are obtained by means of 
a gaussian expansion of integrands. Main integrals are performed with the Simpson's 
method. 

Typical running time: Almost all CPU time is consumed by calculations of integrals, 
specially for transverse momentum distributions which involves multiple integral loops. 
It takes up to 30 min on a 2GHz Intel P4-processor machine. 
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LONG WRITE-UP 



I. INTRODUCTION 

Single-nucleon knockout reactions with heavy ions, at intermediate energies and in in- 
verse kinematics, have become a specific and quantitative tool for studying single-particle 
occupancies and correlation effects in the nuclear shell model The high sensitivity of 
the method has allowed measurements on rare radioactive species available in intensities of 
less than one atom per second for the incident beam. The experiments observe reactions in 
which fast, mass A, projectiles collide peripherally with a light nuclear target, typically ^Be 
or ^^C, producing residues with mass {A — 1), in the following referred to as the core (c) of 
the assumed two-body system of core plus nucleon. The final state of the target and that of 
the struck nucleon are not observed, but instead the energy of the final state of the residue 
can be identified by measuring coincidences with decay gamma-rays emitted in fiight. Re- 
ferred to the center-of-mass system of the projectile, the transferred momentum is k^. In 
the sudden approximation and for the stripping reaction, defined below, this must equal the 
momentum of the struck nucleon before the collision. The measured partial cross sections 
to individual final levels provide spectroscopic factors for the individual angular-momentum 
components j. In complete analogy to the use of angular distributions in transfer reactions, 
the orbital angular momentum / is in the knockout reactions revealed by the distributions 
of the quantity kc. 

The early interest in momentum distributions came from studies of nuclear halo states, 

for which the narrow momentum distributions in a qualitative way revealed the large spatial 

extension of the halo wave function [2j. It was shown by Bertulani and McVoy 3] that the 

longitudinal component of the momentum (taken along the beam or z direction) gave the 

most accurate information on the intrinsic properties of the halo and that it was insensitive 

to details of the collision and the size of the target. In contrast to this, the transverse 

distributions of the core are significantly broadened by diffractive effects and by Coulomb 

scattering. For experiments that observe the nucleon produced in elastic breakup, the trans- 

n 

verse momentum is entirely dominated by diffractive effects, as illustrated :4] by the angular 
distribution of the neutrons from the reaction ^Be(^^Be,^°Be-|-n)X. In this case, the width 
of the transverse momentum distribution reflects essentially the size of the target. 
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The cross section for the production of a given final state of the residue has two con- 
tributions. The most important of the two, commonly referred to as stripping or inelastic 
breakup, represents all events in which the removed nucleon reacts with and excites the 
target from its ground state. The second component, called diffractive or elastic breakup 
[^, represents the dissociation of the nucleon from the residue through their two-body in- 
teractions with the target, each being at most elastically scattered. These events result in 
the removed nucleon being present in the forward beam with essentially the beam velocity, 
and the target remaining in its ground state. These processes lead to different final states, 
they are incoherent, and their cross sections must be added in measurements where only the 
residue is observed. General expressions for the total and_differential cross sections for the 
two components have been given by Hussein and McVoy |6| and further developed by other 
authors in studies of reactions with halo nuclei (see, e.g. ref. 

In a subsequent development, the knockout method was extended to non-halo states 
y, ^, Q, [ill [3, ll^ ll^- For these, involving more deeply-bound nucleons, the one- 
nucleon absorption cross sections are much smaller than the free-nucleon reaction cross 
section on the same target; a ratio that gives a measure of how much the nucleon wave 
function is "shielded" from the target by the bulk of the core. This required a more elaborate 
theo.e«cal .eat.ent b.ed o„ t.e e,a.,c 5-.a.„ce. S. and 5„ E3 of the CO. and 
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nucleon. Other advanced theoretical treatments by Bonaccorso and collaborators 
have shown the relevance of transfer to the continuum in eikonal methods for stripping and 
diffraction dissociation. 

The code described in this article is based on the work presented in ref. j2l| and on 
previous works, as those mentioned above. The code calculates single-particle bound state 
wavefunctions and eikonal S-matrices. Using these wavefunctions and ^-matrices, one ob- 
tains single-particle cross sections after integration over momentum and summation over 
the mi substates. As shown in ref. j^^, the use of this code demonstrates (i) that the 
transverse momentum distributions are quantitatively and even qualitatively different from 
the parallel momentum distributions, and (ii) can serve to extract angular-momentum in- 
formation from the angular distributions of the residues. The code is also of importance 
for calculating acceptance corrections in experiments, and one can evaluate the correlations 



between longitudinal and radial components 
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II. BOUND STATES 



The computer code MOMDIS calculates various quantities of interest for knockout reac- 
tions of the type 

{c + v)+T ^c + X. (1) 

The internal structure of nucleus c, the valence particle v and of the target T is not taken 
into account. The initial state of the projectile nucleus P — c + v is obtained by the solution 
of the Schrodinger equation for the relative motion of c and v in a. nuclear + Coulomb 
potential. Particles c, v, and P have intrinsic spins labeled by Ic, ly and J, respectively. 
The corresponding magnetic substates are labeled by Mc, and M. The orbital angular 
momentum for the relative motion of c + v is described by / and m. In most situations of 
interest, particle v is a, nucleon and c is a "core" nucleus. Thus it is convenient to couple 
angular momenta as 1 + I„= j and j + 1^,= J, where J is the channel spin. Below we also use 
the notation s, instead of I„, for the intrinsic spin of particle v. 
The bound state wavef unctions of P are specified by 

^jM (r) = ^ (r) y'j^ = ^^y'jM , (2) 

where r is the relative coordinate of c and v, uf- (r) is the radial wavefunction and y^j^ is 
the spin-angle wavefunction 

yjM- {3mIcMc\JM)\jm)\IcM,), with |jm) = ^ Xim, (r) Xm. (3) 

m , Mc mi , rris 

where is the spinor wavefunction of particle v and {jmIcMc\JM) is a Clebsch-Gordan 
coefficient. 

The ground-state wavefunction is normalized so that 

oo 

j (fr |^jM(r)|' = jdr \ui^{r)\^ = l. (4) 



The wavefunctions are calculated using a spin-orbit potential of the form 

V{v) = V^{r) + Vs{r) (l.s) + Vc{r) (5) 
where Vo(r) and Vs{r) are the central and spin-orbit interaction, respectively, and Vc{r) is 
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the Coulomb potential of a uniform distribution of charges: 



Vb(r) 



2R, 



c 



for r > Rc 

,2 



3 - 



r 



for r < Rc, 



(6) 



where Zj is the charge number of nucleus i = v,c. 

The nuclear potential is assumed to have a Woods-Saxon form plus a spin-orbit interac- 
tion, 



Voir) = Vo /o(r), and Vsir) 



V< 



so 



h 



Id 
r ar 



with fi{r) 



1 -|- exp 



a,- 



-1 



(7) 



The spin-orbit interaction in Eq. [7| is written in terms of the pion Compton wavelength, 
/i/m^c = 1.414 fm. The parameters Vq, Vso, Rq, o,o, Rso, and aso are adjusted so that the 
ground state energy Eb (or the energy of an excited state) is reproduced. 

The bound-state wavefunctions are calculated by solving the radial Schrodinger equation 

[ d"^ l{l + 1) 



2m, 



ab 



dr'^ 



<• (r) + [Vo (r) + Vc (r) + (s.l) ^50 (r)] ^ (r) = (r) (8) 



where (s.l) = + 1) — /(/ + 1) — s{s + 1)] /2. This equation must satisfy the boundary 
conditions ufj (r = 0) = ufj (r = 00) = which is only possible for discrete energies E 
corresponding to the bound states of the nuclear -t- Coulomb potential. 



III. UNITARITY, STRIPPING AND DIFFRACTION DISSOCIATION 



Using eikonal waves for scattering wave functions, one gets 



^(-)*(r)$(+)(r) = S (6) exp (zq.r) , 



where S (b) is the scattering matrix 

S (6) = exp [ix(^)] 5 with 



X{b) 



1 

hv 



dz Uopt{r) 



(9) 



(10) 



and Uoptij) is the appropriate optical potential for the core+target and the neutron (or 
proton) +target scattering. In equation (fTUj) x(6) is the eikonal phase, and r = Vfe^ + z^, 
where b is often interpreted as the impact parameter. This interpretation arises from a 



comparison of the results obtained with eikonal wavefunctions with those obtained in with 
classical particles colliding with a fixed impact parameter b [2^. Nonetheless, the eikonal 
wavefunction is a quantum scattering state and b is the transverse coordinate associated 
to them. Thus wave-mechanical effects, like smearing and interference, are accounted for 
properly. 

As seen in eq. ^1 in the eikonal approximation, the phase shift is linear in the interaction. 
Then, according to Glauber 2^, the S-matrix for P{c + v) +T can be written as 



Sp — SySc, 



(11) 



where the S is still an operator since it refers to a particle inside the projectile and this 
depends upon parameters that have to be averaged over by the ground state wave function 
of projectile, 0o- The elastic breakup cross section can be easily calculated (the observed 
particle is c) 



O'el.bup — (^k 



SbSc 



where 0k is wave function that represents the continuum of c and v. Since 
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|0k) (0kMk+ |0o) (0o| 
the above expression can be simplified 



Sis: 
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SbSc 


0o) 



(13) 



In so far as the inelastic break up is concerned, one realizes that the detected particle 
must reach the detector intact and thus one must use a survival probability to guarantee 



this. This survival probability is 



Sr 



. On the other hand the interacting fragment v is 

, 2^ 



removed (stripped) and the probability for this to happen is ( 1 

, 2 



Sv 



Identifying the 



Sv 



transmission coefficient with 1- 

(Tin.bup = -j^Yl ('^O 

The above cross-section is also called the v removal cross-section with the notation (T_„. 
Of course one may have the removal of the core (less probable owing to the Coulomb barrier 



we can write for the inelastic break up the following 

(14) 



1 - TJ T„ 



between the charged core and the target) 



TT 



1 - TJ T, 



The sum of cr_^ and cr_c gives 



k 



T(, -\- Ty 2,TyTf, 



(15) 



(16) 



Summing the above with the total absorption of P, described here by cr_p = cr, 



abs 



T T 



0/, gives 



TT V - 



T -\- T —TT 



reaction i 



(17) 



which confirms unitarity. 

The expressions for aei.bup, eq. ESI and (Tin.bup, eq. HH have been used by several authors 
to analyze the data on halo nuclei. It has been common to call aei.bup = tXrfif, the diffractive 



break-up cross section and a.. 



tn.up 



astrip the stripping cross section 
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The cross sections also depend on the single particle occupancy probability which can be 
accounted for multiplying them by spectroscopic factors Si 

IV. EIKONAL S-MATRICES 

In the optical limit of the Glauber theory and the approximation (explained in 

detail in ref. [3|), the eikonal phase is obtained from the nuclear ground state densities and 
the nucleon-nucleon cross sections by the relation (here we drop the operator notation) 

1 



S{b) = exp [ixib)] 



with 



XNib) 



dQ Q Pp{<l) Pt{q) fNN{q) Joiqb) , (18) 



kNN Jo 

where Pp^t (q) is the Fourier transform of the nuclear densities of the projectile and target, and 
Inn (<?) is the high-energy nucleon-nucleon scattering amplitude at forward angles, which 



can be parametrized by 



Inn (g) 



knn 
Air 



aNN (i + cxnn) exp {-^nnQ^) 



(19) 



In this equation cnn, c^nn, and Pnn are parameters which fit the high-energy nucleon- 
nucleon scattering at forward angles. These parameters were originally obtained for Ejsf > 
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100 MeV/nucleon 29|. We use an additional set of parameters which have been shown to 



reproduce the nucleon-nucleon scattering data in the energy region 30 MeV < En < 2200 
MeV 3]. The nucleon-nucleon amplitude fNN{(l) which enters equation [THl is of course 
not the same as the one expressed in equation for the free nucleon-nucleon scattering. 
Medium modifications (e.g., isospin average, etc.) have to be taken into account. As shown 



in ref. j28|, the largest medium effect is due to Pauli-blocking |30]. The code MOMDIS 
allows (optional) one to include Pauli blocking effects in the nucleon-nucleon cross section 
a^N, according to refs. Q] and 3^ . 

In eq. (fT^ the quantities pp (g) and pt (g) are calculated from the radial density distri- 
butions, usually taken to be of Gaussian shapes for light nuclei, and of Fermi shapes for 
heavier nuclei with parameters taken from experiment. For cases where more accuracy is 
needed, it is possible to take the density distributions directly from Hartree-Fock calcula- 
tions. A set of experimental values of density parameters was published in ref. . These 
are useful, in particular in parametrizing the densities of stable nuclei. 

For the Coulomb part of the optical potential the integral in eq. (fTIHl diverges. One solves 
this by using x = Xn + Xc, where xn is given by eq. (finjl without the Coulomb potential 
and writing the Coulomb eikonal phase, xc as 

Xc{b) = 2r]\n{kb) , (20) 

where rj = ZpZte'^/hv, Zp and Zt are the charges of projectile and target, respectively, v is 
their relative velocity, k their wavenumber in the center of mass system. Eq. (j^UI) reproduces 
the exact Coulomb scattering amplitude when used in the calculation of the elastic scattering 
with the eikonal approximation |22l |: 

fc{0) = , f'^tl/,, {-^V In [sm\e/2)] + .vr + 2^0o} (21) 
2pv^ sm [6/2) 

where 0o = argr{l + iri/2). This is convenient for the numerical calculations since, as shown 
below, the elastic scattering amplitude can be written with the separate contribution of the 
Coulomb scattering amplitude included. Then, the remaining integral (the second term on 
the r.h.s. of eq. (j26|) below) converges rapidly for the scattering at forward angles. 

Although the Coulomb phase in eq. ((201) diverges at 6 = 0, this does not pose a real 
problem, since the strong absorption suppresses the scattering at small impact parameters. 
It is also easy to correct this expression to account for the finite charge distribution of the 



nucleus. For example, assuming a uniform charge distribution with radius R the Coulomb 
phase becomes 



Xc{b) = 27] \ 6(6 - R) \n{kb) + e{R - b) \n{kR) + ln(l + VI - b^R^) 

1 



'1-byR^- -{l-b'/R 
3 



2\3/2 



(22) 



where B is the step function. This expression is finite for 6 = 0, contrary to eq. EB If one 
assumes a Gaussian distribution of charge with radius R, appropriate for light nuclei, the 
Coulomb phase becomes 



Xcib) = 2r] 
where the error function Ei is defined as 

Ei{x) 



Hkb) + -E,{byR') 



(23) 



oo g-t 



dt . 



(24) 



This phase also converges as 6 — * 0. In eq. R = Rp + Rt, while in eq. R 



Rp + Rf, where Rp and Rt are the respective projectile and target radius. The cost of 
using the expressions (j^^ and is that the Coulomb scattering amplitude becomes more 
complicated than ()21|) . Moreover, it has been proven 22] that the elastic and inelastic 
scattering cross sections change very little by using eqs. or (j2ni), instead of eq. (j^ . 

In calculations involving stripping the final state Coulomb interaction between the core 
and the target is taken into account by using the eikonal-Coulomb phase shift of ()2()|1 in the 
calculation of Sc- However in the calculation of diffraction dissociation both Sc and Sn are 
calculated using the eikonal-Coulomb phase shift of 



V. ELASTIC CROSS SECTIONS 



The calculation of elastic scattering amplitudes using eikonal wave functions, eq. 



is 



very simple. They are given by [22 

fel{0)=tk 



db b Joiqb) {1 — exp [«x(^)]} 



(25) 



where q = 2ksin{6/2), and 6 is the scattering angle. The elastic scattering cross section is 
daei/dQ = \ fei{0)f. For numerical purposes, it is convenient to make use of the analytical 
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formula, eq. (EH), for the Coufomb scattering amplitude. Thus, if one adds and subtracts 
the Coulomb amplitude, fc{d) in eq. (f^ . one gets 

/•oo 

felie) = fc{0) +ik dbb Mqb) exp [txcib]] {1 - exp . (26) 

Jo 

The advantage in using this formula is that the term 1 — exp [iXN{b)] becomes zero for 
impact parameters larger than the sum of the nuclear radii (grazing impact parameter). 
Thus, the integral needs to be performed only within a small range. In this formula, xc is 
given by eq. (plj) and fc{0) is given by eq. (PT|). with 

00 = -VC + Y] I - arctan —— ] , (27) 

and C = 0.5772156... is the Euler's constant. 

The elastic cross section can be expressed in terms of the transverse momentum by using 
the relationships dQ ~ d'^k±/k'^, and k± c:^ q = 2k sm{6/2), valid for high-energy collisions. 



VI. MOMENTUM DISTRIBUTIONS 



Following ref. j6[ the stripping reaction {c + v) + A — > c + X, where c corresponds to a 
specified final state of the core, is given by 

2 



da. 



str 



d'kc (27r) . - ^ 
and where r = (p, z, 0) = — r^, so that 

bc = \p- b,„| 



(28) 



+ bl-2p b^ cos I 



sin^ 9 + bf, — 2r sin 9 b^ cos I 



(29) 



Assuming an independence of the bound state wave functions on spin (i.e. no spin-orbit 
term), equation ESI can be recast to 

2 



r e 



-ikc.r 



Sc (be) fpim (r) 



(30) 



d^k, (27r)= 

where the single-particle bound state wave functions for the subsystem (c + f ) is specified by 
i'lm (i") = Ri (^) (r), where Ri (r) is the radial wave function. The code MOMDIS uses 
eq. Eni (and the subsequent ones below) to obtain the momentum distributions. But it also 
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allows one to calculate and use radial wavefunctions with a spin-orbit potential term which 
is useful in some cases. 

The cross sections for the longitudinal momentum distributions are obtained by integrat- 
ing eq. (j28|) over the transverse component of kc, i.e. over k^, and using 



J d'kj exp [-tk^p-p]] = {2nf6{p~p') . 



(31) 



One gets 



da. 



str 



dh (27r) 



1 1 roo poo 

^^T^J^l [i-\SAb.)f] d'p\s^ 



(Mr 



X 



dz exp [-ik^z] ipim (r) 



(32) 



where represents the longitudinal component of k^. 

For the transverse momentum distribution in cylindrical coordinates k_i_ = ^Jk'^ + k'y, one 
uses in eq. f^H]) 

/oo 
dkz exp [—ikz{z — z')] = 2n5 {z — z') , (33) 
-oo 

and the result is 



da. 



str 



d'^ki, 27r2/ + l Jq 



d\, [l-\S,{K)f] 



/OO n 
dz I d^p exp [-ikj.p) (K) ipim [r] 



(34) 



The code MOMDIS also allows to calculate the transverse momentum distributions in 
terms of the projection onto one of the Cartesian components of the transverse momentum. 
This can be obtained directly from eq. (jHU, i.e. 

da^tr f dastr 



dkri, 



dkr 



d'^ki 



{kxi ky) . 



The total stripping cross section can be obtained by integrating either eq. 
fj3l|l . For example, from eq. (j32|l . using eq. (j33|) . one obtains 

27r '"^ 



(35) 
or eq. 



C"str 



21 



poo 

— / db^ by [l - \Sy (6i,)|^] 
^ Jo 



X / d^^r 



Sc 1 \/ sin 6 + bl — 2r sin 6 6„ cos i 



(36) 
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Using the explicit form of the spherical harmonics 



mj2l + ll il-m)\ 



= Ci^Pi^ {cos 6) e'^^ (37) 

and 

kc • r = kx^r sin 6 cos {(pk — 0) + k\f cos 6 , (38) 
part of the integral in ()28|) is shown to be of the form 

J'lrn {k^, k,, K) = [ d\ e-*^-^^e (&c) (r) 



Cim J dr sin 9 dO dcf) exp | —i [A;_|_r sin 6 cos (</)fc — 0) + k\\r cos 6'] } 

(39) 



X Sc \/r2 sin^ ^ + 52 _ 2r sin 9 K cos (0 - 0„) LR, (r) P^^ (cos ^) e^™"^ 



To simplify the calculations we can express Sc (&c) as an expansion in terms of integrable 
functions. The S-matrices can be well described by the expansion 

Sc (be) = Yl ' ^ith = — ■ (40) 

Good fits for realistic S-matrices are obtained with = 20, i.e. with 20 complex coefficients 
aj and Rl = 10 — 20 fm, depending on the size of the system. Since the real part of the 
S-matrices has the usual behavior of Sc (be) ~ for be ^ R, and Sc (be) ~ 1 for be ^ R, 
where i? is a generic nuclear size, one of the coefficients of the expansion in eq. ()40|) is 
aj = 1, and (3j = 00, which we take as the j = term in the expansion. 

The expansion EUl allows to calculate analytically many of the integrals in eq. ()28|) . and 
one gets (see Appendix of ref. |l21l] ) 

k±dkj_dkz 21 + 

where 



db^ b^ [1-15^ {bv)f] ^ Cf^ \Aimp {kx, k-,, b^ 

m, p 



(41) 



Aimp{k±,k^,bv) = exp [-bl/P]] 

dp p Jp {kj_p) exp [-p'^/P]] I, 



1 

CO 



m—p 



2pbv 
PI 

/oo 
dz exp[-iA;^z] Ri{r) Pi^n (cos 6) . (42) 
-00 
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The first term of tlie Equation (jl^ . witli Pj = oo and aj = 1 can be calculated using 

la (0) = 5.. 

Using the integral of eq. in eq. (PT|) one gets for the transverse momentum distribu- 
tion 



where 



—— / dKK Y.^L / dz \Vi^,{k^X.z)\\ (43) 



N 



X y (ip p Jp (A;_Lp)exp [-pV/?j] -^m-p J -R/ (?^) Pim{cos9) . 

The momentum distributions from diffractive dissociation have nearly the same shape as 



those from stripping j2l|. The code MOMDIS uses the equations of this section to obtain 
the momentum distributions renormalized so that their integral over momentum yields the 
total knockout cross section, eqs. ITSlplusIT^ 

For cases where the magnetic quantum number m differs from zero, we weight the dif- 
ferential cross section with the multiplicity of 2, so that the sum over all m components 
gives the total cross section. The two-dimensional momentum distribution does not depend 
on the azimuthal angle. It is convenient to present it as a function of the parallel and the 
transverse momentum with the definition 

dk±dkz d'^k±dkz^ 

which normalizes to the total cross section when the integration is extended over the negative 
k^ axis. 



VII. COMPUTER PROGRAM AND USER'S MANUAL 
A. Distribution 

The package is distributed in a tar.gz file and, under UNIX systems, can be unpacked as 
follows: 

gunzip filenameXax .gz 
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tar -xvf filenameXax 

Unpacking the file generates two files, a README file, the Fortran source, 
MOMDIS.FOR, and a directory, testdata. This directory contains sample input files for 
a typical run of MOMDIS and files containing the expected output. 

Option 1 - Eigenfunctions and energies 

input file, a_bound.txt - output file, bound. out 

Option 2 - S matrices - Sn 

input file, a_sn.txt - output file, s_sn.out. 

Option 2 - S matrices - Sc 

input file, a_sc.txt - output file, s_sc.out 

Option 3 - Momentum distributions 

input file, a_bound.txt, bound. out, s_sn.out, s_sc.out 

output file, sigma.out 

The MOMDIS.FOR program file is in DOS format. To use on a unix system it should 
first be converted. 

dos2unix -ascii MOMDIS.FOR MOMDIS.f 

It can be compiled using the command 
g77 -o MOMDIS.exe MOMDIS.f 

and run using the command 
./MOMDIS.exe 

The screen output produced from each of the sample runs mentioned above is shown in 
the README file contained in the distribution. 

B. Fortran code 

The units used in the program are fm (femtometer) for distances and MeV for energies. 
The output cross sections are given in miUibarns (mb), mb/(MeV/c), and mb/(MeV/c)^. 

The program is very fast, except for transverse and double-differential momentum distri- 
butions, and does not require a complicated input. It is divided in 3 modules: 
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Module 1 - Calculates of energy and wavefunction of bound states. Module 2 - 
Calculates S-matrices. Module 3 - Calculates momentum distributions. If one chooses 
option Module 3, other options follow: 1 - for longitudinal momentum distributions, eq. 1221 
2 - for transverse momentum distributions, eq. 1321 3 - for projected transverse momentum 
distributions, eq. ESI 4 - for double differential momentum distributions, eq. ESI and 5 - 
for elastic scattering cross section of the core with the target, using eq. EHl 

C. Sample input file 

A sample input file is shown below. It is prepared to study the reaction ^^C+ ^Be — *• 
^^C +X at 103 MeV/nucleon laboratory energy. The rows starting with a symbol "*" are 
not read as input and can be used in the input file at free will. They are well suited to remind 
the user of the input procedure. For more details on how to prepare an input file see the 
readme file located at the end of the Fortran code. With the sample input below the code 
calculates the bound-state wave function. To obtain S-matrices, one needs to comment with 
a "*" the input lines (4th and 5th rows) for bound-state wavefunction and delete the "*" 
from the input lines for S-matrices. The comments in the following two input lines depend 
upon the option parameters Ismat , ^dproj, and Idtarg- For momentum distributions, 
only the first three rows of the sample file below are used. 
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A more specific description follows. Row 1: Zp, Ap = projectile charge and mass 
numbers, Zt, At = target charge and mass numbers, E^ab = laboratory energy per nucleon 
(MeV). Row 2: Zcore, ^core = projectile core (^^C) charge and mass number. Row 3: 
S^C = Spectroscopic factor. Row 4: Single-particle bound state quantum numbers, hq 
= nodes of the wave function (exclude origin), jo = total angular momentum, Iq = orbital 
angular momentum. Row 5: Parameters of the Woods-Saxon potential used to calculate 
wavefunction following eq. El Vq = depth of central potential, Vso = depth of spin-orbit 
potential, Rq = radius parameter of the central potential, A = diffuseness of the central 
potential, Rso = radius parameter of the spin-orbit potential, Aso = diffuseness of the spin- 
orbit potential, Rc = Coulomb radius parameter (usually, Rc = Rq)- Row 6: Ismat = 0, 
or 1. If Ismat = 0, S-matrix is for valence particle-|-target. If Ismat = 1, S-matrix is for 
core+target. Row 7: IpoT = 0, 1, or 2. If IpoT = 0, optical potential is a Woods-Saxon, 
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with real and imaginary parts, V + iVj. If IpoT = 1? optical potential is entered from an 
external, user provided, file. If IpoT = 2, the optical potential is built with the t-pp method. 
Row 8: If IpoT = 0, enter Vq [Vj] = real part [imaginary] (both > 0) of Woods-Saxon 
potential. Rq [Roi] = radius parameter {R = Rq{A](^ + AI/^). d [di] = diffuseness. Row 
9: Enter names of the input files for the real (FILEl) and imaginary (FILE2) parts of the 
potential. In the input files, the first row must give the number of the remaining rows which 
form an ordered list of r (in fm) x V [Vj] (in MeV) in steps of a constant value of Ar. 
Row 10: If IpoT = 2, enter the following options: IpAULi = for no Pauli correction of nn 
cross sections, Jpyi[7i7- = 1 for Pauli correction of NN cross sections. Row 11: Enter options 
for projectile and target densities. Idproj, Idtarg = 1,2,3,4,5. If Idproj = 1; projectile 
density is a three-parameter Gaussian density 

,2 



G{r) = l + b(^) exp ( ) . (46) 



If Idproj = 2, projectile density is a two-parameter Yukawa 



r(r) =r'^-^exp(-/5r). (47) 



If Idproj = 3, projectile density is a four-parameter Fermi function 

1 



F r) 



1 + M- 

a 



(48) 



If Idproj = 4, projectile density is a power of the Fermi function {Ider = 0), or its 
derivative {Ider = 1) 

S{r) = —. --^ , or — . (49) 

[l + exp(^)]" ' dr 

If Idproj = 5, projectile density is a calculated from liquid-drop model. If Idproj = 10, 

projectile density is from an input file. Row 12: In case Idproj = 1 or 2, enter here 

Gaussian or Yukawa density parameters for the projectile, according to eqs. HHl and W7\ 

respectively. Row 13: If Idproj = 3 enter the four parameter of the Fermi function, eq. 

EHl Row 14: If Idproj = 4, enter Fermi density parameters for the projectile and its power, 

as in eq. EHl Row 15: If Idproj = 10, enter name (FILES) of the input file where the 

projectile density is to be found and ap is the proton size parameter to be used in a folding 

of the input distribution with the proton gaussian density, Pp{r) = (vra^) ^''^ exp (— r^/a^). 

In the input file FILES, the first row must give the number of the remaining rows which 
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form an ordered list of r (in fm) x ppROj (in fm~^) in steps of a constant value of Ar. 
Rows 16, 17, 18 and 19: Same as the previous four rows, but for the target (depending 
upon the option Idtarg)- 

D. Eigenfunctions and energies 

Option 1 calls the subroutine EIGEN. The potential parameters to build the Woods- 
Saxon are entered in the input file. In the sample file above, the ground state of ^^C is 
calculated. It is assumed to be a lsl/2+ s-wave {n — 1, I — j — 0.5). This is just one of 
the possible single particle occupancies in the ground state of ^^C. The spectroscopic factor 
can be adjusted to account for the occupancy probability of this single-particle configuration. 

The calculations are mainly done in the subroutine BOUNDWAVE which solves the 
Schrodinger equation for the bound-state problem. When Woods-Saxon potentials are used 
they are constructed in the routine POTENTIAL. 

The output of the wavefunction will be printed in EIGEN.TXT and another file specified 
by the user, which will be ready for use later. 

E. S-matrices 

Option 2 calls subroutine SMAT. If this option is used the input rows for the bound- 
state wavefunction should be commented with a "*" . Here one has to specify if one wants 
Sy or Sc and if the optical potentials are entered in a separate file or built with the "t-pp" 
approximation. If the later is chosen, one has to enter the options for the nuclear densities. 
An option is given to include, or not, the effect of Pauli blocking on the in- medium nucleon- 
nucleon cross section. 

F. Momentum distributions 

Option 3 calls for calculation of momentum distributions. In this case, only the general 
input (charges, masses and bombarding energy) are used. All subsequent rows may be 
commented, or not. One can then choose to calculate da/dpz, (Pa/dp^, da/dpy, d^a/dptdpz, 

or CTelast- 
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VIII. ROUTINES INCLUDED WITH THE CODE DISTRIBUTION 



EIGEN finds eigenvalue and eigenfunction of a system composed of a valence particle 
with charge Zl and a nucleus with charge Z2. The nuclear potential is a Woods-Saxon, with 
spin-orbit and Coulomb components. 

BOUND WAVE solves the radial Schrodinger equation for bound states. 

POTENT constructs potentials to be used in the calculation of bound state wavefunc- 
tions. Uses routine OMP_WS which constructs the Woods-Saxon potential. 

SMAT computes the eikonal S-matrices for heavy ion collisions. The S-matrices are 
functions of the impact parameter b. The program is appropriate for heavy ion bombarding 
energies EJab > 30 MeV/nucleon. 

CROSS. If lOPT = 1, this routine calculates da /dp,. If lOPT = 2, 3 , 4 and 5, it 
calculates d'^a/dp'^, da/dpy, d'^a/dptdpz, or (Teiasti^), respectively. 

COMPDIF is used in the calculation of the total diffraction dissociation cross section, 
eq. d 

GAUST calculates the sum over the Gaussian expansion for da/dpy, where Py is one of 
the components of the transverse momentum . 

GAUSLT is similar to GAUST, but calculates d^a / dptdpz- 

FINTEG calculates dz \Ri (r)|^ Pi^ {cos 6) which is used in several other integrals. 
PINTEG deals with the integral j d(p S (^^/r2 sin^ 6 + - 2r smO h cos 0^ ^ . ZINTEG 
computes dz exp [—ik^z] Ri (r) Pim (cos 6*). COMPL is used in the calculation of the 
longitudinal momentum distribution, eq. ESI 

SMAT_FIT finds the best expansion parameters aj to fit Sc with an expansion in gaus- 
sians (eq. Uses auxihary routines FUNGS, LFIT, GAUSS J and COVSRT. 

FCOUL and CPO are used in the calculation of the Coulomb elastic amplitude, eq. |^ 

SIGNNE and PHNNE are used in the calculation of the nucleon-nucleon scattering 
amplitude, eq. UHl 

PHNUC calculates the eikonal nuclear phase, eq. UHl 

The routines GAUSS (eq. HHl), YUKAWA (eq. ED, FERMI (eq. HHD, SAXON 
(eq. and DROP are used to calculate the nuclear densities parametrized as Gaussians, 
Yukawa, Fermi, Woods-Saxon, and liquid drop functions, respectively. The parametriza- 
tions follow the equations 1^611^ DROP builds liquid drop model densities for the nuclei. 
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following ref . 32 1 . 



The routines PLGNDR, BESSJ, BESSJO, BESSJl, BESSI, BESSIO and BESSIl 

calculate Legendre polynomials and Bessel functions. 

The other routines in the code are shortly described as follows. GFV calculates factori- 
als, BETWEEN and FINFOUT are interpolating routines, SKIPCOM is used to skip 
comments in the input file, RSIMP does Simpson's integrations of tabulated functions, 
FOURIERO calculates Fourier transforms (e.g. in the calculation of the transformed den- 
sities of eq. ITHj) . NUCNAME finds the symbol associated to a given chemical element and 
EXCESS calculates the mass excess with the parameters published in ref. 

IX. THINGS TO DO 



1 - Use the sample input file to obtain fig. 6 of ref. |21| 

2 - Modify the input file so as to reproduce all other results published in ref. |21| |. 

3 - Use the overlap function for "^jm (i") published in ref. 2^ to reproduce their figure 2. 
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